* uses the following packages from SSC: coefplot, estout, blindschemes, reghdfejl, julia
* if "reghdfejl" calls are replaced with "reghdfe" then reghdfejl and julia are not needed but reghdfe is

cap cd "/Users/davidroodman/Library/CloudStorage/OneDrive-OpenPhilanthropyProject/Replication opinions/Heyes and Saberian 2019"
cap cd "D:/OneDrive - Open Philanthropy Project/Replication opinions/Heyes and Saberian 2019"

set scheme plottig
set varabbrev off

if c(os)=="Windows" global font LM Roman 9
               else global font Latin Modern Roman
graph set window fontface "$font"


***
*** HS19 & HS22 regressions
***

global weather6t4 temp6t4 press6t4 dew6t4 prcp6t4 wind6t4 skycover
global pollutants ozone co pm25
global dummies dow nat_name c_asy_type year ct#month chair
forvalues i=1/2 {
  local yr  : word `i' of 19 22
  local if  : word `i' of 1 `"inlist(c_asy_type,"E","I") & gender!="""'
  use "Heyes and Saberian 20`yr'/Data/final/matched" if `if', clear
  cap rename pm pm25
  gen year = year(date)
  gen month = month(date)
  gen dow = dow(date)
  gen m = mofd(date)
  encode city, g(ct)
  eststo HS`yr'_0004CH  : reghdfejl res $weather6t4 $pollutants if nat_name=="CHINA", cluster(ct#month) a($dummies) keepsing
  eststo HS`yr'_0004NoCH: reghdfejl res $weather6t4 $pollutants if nat_name!="CHINA", cluster(ct#month) a($dummies) keepsing
  eststo HS`yr'_0004    : reghdfejl res $weather6t4 $pollutants                     , cluster(ct#month) a($dummies) keepsing
}

* HS22 preferred estimate but without pollution controls, cited in opinion
reghdfejl res $weather6t4 if nat_name!="CHINA", cluster(ct#month) a($dummies) keepsing

* inspired by Fischman (2014), Figure 8--showing cross-judge variation in NYC for Chinese
preserve
gen Chinese = nat_name=="CHINA"
collapse res (count) N=res if hearing_loc_code=="NYC", by(chair Chinese)
keep if N>=50  // per Fischman, only keep means with N>=50
drop N
reshape wide res, i(chair) j(Chinese)
scatter res1 res0, xsize(5) ysize(5) ytitle(Chinese applicants) xtitle(non-Chinese applicants, margin(small)) ///
  graphregion(margin(0 4 0 2)) plotregion(margin(zero)) name(grantrates, replace) ///
  xlabel(0 "0%" .2 "20%" .4 "40%" .6 "60%" .8 "80%" 1 "100%") ylabel(0 "0%" .2 "20%" .4 "40%" .6 "60%" .8 "80%" 1 "100%") ///
  note("Notes: Each dot represents a judge. Copying Fischman (2014, FIgure 8), only New York judges with" ///
       "at least 50 decisions are included. The HS22 data span January 1, 2000—September 30, 2004.", size(vsmall) margin(medsmall) xoffset(-11))
graph export replication/grantrates.png, width(2000) replace
restore


***
*** S22 regressions
***

global weather6t4 temp6t4 prcp6t4 skycover press6t4 dew6t4 wind6t4
global noaaweather temp_noaa_6a4p prec_noaa_6a4p sky_noaa_6a4p pressure_6a4p dew_noaa_6a4p wind_speed_noaa_6a4p
global pollutants pm25 ozone carbon_monoxide
global dummies judge_B defensive nationality dow fyear ct#month \$mid
global mid

cap frame drop weather
mkf weather
frame weather: use location_code date $noaaweather $pollutants using "Spamann 2022/data/weather/NOAA_AQS" if !mi(temp_noaa_6a4p) & substr(location_code,5,.)=="base"
frame weather: replace location_code = substr(location_code,1,3)
frame weather: rename ($noaaweather) ($weather6t4)

use base_city_B comp_date_B latest_hearing_A grant judge_B defensive nationality hearing_loc_B using "Spamann 2022/data/EOIR_asylum/asylum_revised" if !mi(comp_date_B), clear // cf. asylum_data_cleaning_and_checking.do: already limited to years 1990-
decode base_city_B, gen(location_code)
encode location_code, gen(ct)
gen int year = year(comp_date_B)  // note in S22, these always depend on comp_date, even when using last_hearing elsewhere
gen int fyear = year + (month(comp_date_B)>9)

gen int date = latest_hearing_A if latest_hearing_A<=comp_date_B
frlink m:1 location_code date, frame(weather)
frget $weather6t4 $pollutants, from(weather)
gen byte month = month(date)
gen byte dow = dow(date)
egen city_month = group(ct month)

* S22 preferred but with two-way clustering & without pollution variables and missingness issue, cited in opinion
reghdfejl grant $weather6t4 if nationality!="CH":nationality & inrange(comp_date_B,td(1jan2000),td(30sep2004)), a(city_month judge_B defensive nationality dow fyear) cluster(ct) keepsing


* generate missingness dummies
foreach var in $pollutants {
  gen byte mid_`var' = mi(`var')
  global mid $mid mid_`var'
  recode `var' (.=0)
}

eststo S22_0004    : reghdfejl grant $weather6t4 $pollutants if  inrange(comp_date_B,td(1jan2000),td(30sep2004))                                , a($dummies) cluster(ct) keepsing
eststo S22_0004NoCH: reghdfejl grant $weather6t4 $pollutants if  nationality!="CH":nationality & inrange(comp_date_B,td(1jan2000),td(30sep2004)), a($dummies) cluster(ct) keepsing
eststo S22_0004CH  : reghdfejl grant $weather6t4 $pollutants if  nationality=="CH":nationality & inrange(comp_date_B,td(1jan2000),td(30sep2004)), a($dummies) cluster(ct) keepsing
eststo S22_9019    : reghdfejl grant $weather6t4 $pollutants                                                                                    , a($dummies) cluster(ct) keepsing
eststo S22_9019NoCH: reghdfejl grant $weather6t4 $pollutants if !(nationality=="CH":nationality & inrange(date,td(1oct1996),td(10may2005)))     , a($dummies) cluster(ct) keepsing  // CPC was enacted 1996/9/30 (PL 104-208) and the conditionality ended 2005/5/11 (PL 109-13)
eststo S22_9019CH  : reghdfejl grant $weather6t4 $pollutants if   nationality=="CH":nationality                                                 , a($dummies) cluster(ct) keepsing

* Like S22_9019NoCH, but excluding *all* Chinese, as in S22, Table 1, col. 10. Doesn't quite match S22. Mentioned in opinion footnote.
reghdfejl grant $weather6t4 $pollutants if nationality!="CH":nationality, a($dummies) cluster(ct) keepsing

coefplot HS19_0004 S22_0004 HS22_0004, bylabel(2000–04) || ///
          HS19_0004NoCH S22_0004NoCH HS22_0004NoCH, bylabel(2000–04, excluding Chinese) || ///
          HS19_0004CH S22_0004CH HS22_0004CH, bylabel(2000–04, Chinese only) || ///
          (,drop(*)) S22_9019 (,drop(*)), bylabel(1990–2019) || ///
          (,drop(*)) S22_9019NoCH (,drop(*)), bylabel(1990–2019, excluding Chinese) || ///
          (,drop(*)) S22_9019CH (,drop(*)), bylabel(1990–2019, Chinese only) ///
          keep(temp6t4) rescale(1000) cismooth(n(10) lwidth(3 3)) ///
          ylab(.75 "Original (HS19)" 1 "Comment (S22)" 1.25 "Rejoinder (HS22)", notick labcolor(black) labsize(medsmall)) ///
          xtitle(Impact on probability of asylum grant (percentage points per 10°F)) ///
           mlabel(string(@b,"%9.3f")+" ("+string(@se*1000,"%9.2f")+")") mlabpos(12) mlabcolor(black) mlabsize(medsmall) ///
          byopts(noiytick legend(off) graphregion(margin(zero)) cols(3) ///
            note("Notes: Point estimates and 95% confidence intervals shown. Standard errors in parentheses, clustered by city-month (HS19 and HS22) or city (S22). Most" ///
                 "displayed results do not appear exactly in HS19, S22, or HS22, but are obtained by modifying public code to restrict the sample as indicated, or, for S22," ///
                 "to cluster in that paper's preferred way. Chinese applications excluded only for October 1, 1996–May 10, 2005.", size(vsmall) margin(medsmall))) ///
          name(asylum, replace)
graph export replication/asylum.png, width(2000) replace


***
*** parole
***

use "Heyes and Saberian 2022/Data/final/parole", clear 
eststo parole1215: reghdfejl res tempmean10 rhmean prcp pressmean windmean ozonemean comean nomean, cluster(insm) a(insm year dayofweek other chair type) keepsing

use "replication/parole1619", clear 
eststo parole1619: reghdfejl res tempmean10 rhmean prcp pressmean windmean ozonemean comean nomean, cluster(insm) a(insm year dayofweek other chair type) keepsing

mat b = -.974  // HS19, Table 7, col. 1
mat se = .482
mat colnames b = tempmean10

coefplot (mat(b), se(se)) parole1215 (, drop(*)), bylabel(2012–15) || (, drop(*)) (, drop(*)) parole1619, bylabel(2016–19) keep(tempmean10) cismooth(n(10) lwidth(3 3)) ///
          ylab(.75 "Original (HS19)" 1 "Rejoinder (HS22)" 1.25 "This opinion", notick labcolor(black) labsize(medsmall)) xtitle(Impact on probability of parole grant (percentage points per 10°F)) ///
          byopts(noiytick legend(off) note("Notes: Point estimates and 95% confidence intervals shown. Standard errors clustered by venue-month in parentheses.", size(vsmall) margin(medsmall)) graphregion(margin(zero))) ///
          mlabel(string(@b,"%9.3f")+" ("+string(@se,"%9.2f")+")") mlabpos(12) mlabcolor(black) mlabsize(medsmall) name(parole, replace)
graph export replication/parole.png, width(2000) replace
